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We present the results of a computer simulation study of thermodynamical properties of TIP4P 
water confined in a hydrophobic disordered matrix of soft spheres upon supercooling. The hydrogen 
bond network of water appears preserved in this hydrophobic confinement. Nonetheless a reduction 
in the average number of hydrogen bonds due to the geometrical constraints is observed. The liquid 
branch of the spinodal line is calculated from 350 K down to 210 K. The same thermodynamic 
scenario of the bulk is found: the spinodal curve is monotonically decreasing. The line of maximum 
density bends avoiding a crossing of the spinodal. There is however a shift both of the line of 
maximum density and of the spinodal toward higher pressures and lower temperatures with respect 
to bulk. 
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I. INTRODUCTION 

The behaviour of the liquid-gas spinodal in water plays 
an important role in the interpretation and the determi- 
nation of its complete phase diagram. In particular the 
related behaviours of the hquid branch of the spinodal 
and the line of maximum density (TMD line) in the su- 
percooled liquid state are connected to the possible exis- 
tence of a second critical pointiiS,. Among the different 
interpretations of the anomalies of water at low temper- 
atures in fact a number of experimental and computer 
simulation results and theoretical calculations indicate 
the possibility of a liquid-liquid (LL) coexistence and a 
second critical point^i^. The low density and high density 
liquid phases would be in correspondence with the low 
density and high density amorphous forms of water be- 
low its glass transition temperature. In accordance with 
these hypotheses the many numerical studies performed 
to determine the spinodal line of bulk supercooled water 
agree in finding that this line is monotonically decreasing 
at least down to the lowest temperatures that was pos- 
sible to investigate^iii^i^iii^'^ii^'i^. Correspondingly the 
TMD line bends to avoid crossing the spinodal. This sce- 
nario excludes both TMD and spinodal behaviour to be 
responsible of the water singularities. These singularities 
can therefore be connected to the existence of the second 
critical point or to a singularity free scenarioi^. 

An increasing amount of theoretical and experi- 
mental studies have been performed in recent years 
on water in confined geometries, see for exam- 

p^ glS. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32 Water is 

in fact confined in many situations of interest for bio- 
logical and technological applications. In this respect 
the study of the modification of the phase diagram of 
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water induced by the perturbation of the confining envi- 
ronments is particularly significant. The study of water 
in confinement also can help to shed light on the be- 
haviour of the bulk in the region of deep supercooling 
were nucleation prevents experiments to be performed 
in the bulk while permits very low temperatures to be 
reached in confinemen t^^ '^"^i^^ displaying scenarios that 
can be connected to the bulk behaviour—. The shift of 
the phase diagram due to confinement could drive the 
so far experimentally inaccessible zone where the second 
critical point is supposedly located in a region accessi- 
ble to experiment. A study of ST2 water confined be- 
tween smooth plates has indicated the possibility of a LL 
transition in confinement^. A recent theoretical work 
indicates a shift to lower temperatures, higher densities 
and higher pressures of the second critical point for hy- 
drophobic confinement^^. Indications of a LL transition 
have been also found in a recent computer simulation 
studjs22, on TIP5P water confined between hydrophobic 
plates, a shift to lower temperature of the second critical 
point is also estimated. 

Generally speaking solid interfaces strongly distort the 
HB network^^. In some cases an hydrophobic environ- 
ment for water can be provided by the presence of an 
apolar solute, as in the case of aqueous solutions of rare 
gases, polymers or large organic molecules. When the 
hydrophobic units are small enough it is expected that 
water maintains its hydrogen bond structure^^^ in spite 
of the perturbation of the solute, this leads to the forma- 
tion of water cages around the non polar solute. In this 
paper we explore how this hydrophobic effect influences 
the behaviour of the spinodal line and the TMD upon 
cooling in the limit of small hydrophobic units to study 
the case where the structure of water is preserved or only 
slightly modified. Since the solute-solute interaction is 
expected not to play an important role in comparison 
with the hydrogen bond and the solute-solvent interac- 
tions, in the present study we consider as a first approxi- 
mation a model where few solute particles are kept fixed 
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to form an hydrophobic confining matrix for water. We 
expect this model to catch the main features of the phe- 
nomenon. 

II. MODEL AND COMPUTER SIMULATION 
DETAILS 

The water fluid is simulated with the DL_POLY pack- 
age (W. Smith, T. R. Forester and I. T. Todorov, Dares- 
bury Laboratory, UK) using the TIP4P site model'^^. 
In this model three sites are arranged according to the 
molecular geometry. The two sites representing the hy- 
drogens are positively charged with qn = 0.52 atomic 
charges, each one forms a rigid bond with the site of 
the oxygen at distance 0.9752 A. The angle between the 
bonds is 104.52°. The site of the oxygen is neutral while 
a fourth site carries the negative charge of the oxygen 
qo = —2qH- This site is located in the same plane of 
the molecule at distance 0.15 A from the oxygen with an 
angle 52.26° from the OH bond. The intermolecular in- 
teractions are represented by Coulombian terms between 
the charged sites and a Lennard- Jones (LJ) potential be- 
tween the neutral oxygen sites. The LJ parameters are 
given by ao = 3.16 A and eo — 0.65 kJ/mol. 

Water is embedded in a rigid disordered array of Nm = 
6 soft spheres. The soft spheres interact with the water 
oxygen sites by means of the potential 

voM{r) = AeoM (~^) (1) 

where 

eoM = -s/eoEM (2) 

and 

com — -{co + ^m) (3) 

We assume ctm = '2<to and eM — O.leo- The average 
distance between the centers of two soft spheres is 10.5 A. 
With this size the soft spheres are far from the limit where 
they appear as hard wall to water. On the other end 
to enhance the volume excluded effects we use a pure 
repulsive water-spheres interaction where the parameter 
Em of the soft potential is taken as O.leo to soften the 
repulsive ramp. The order of magnitude of the size of 
the spheres is close to the typical hard sphere diameter 
of apolar solutes^i*^. 

The simulations are performed in the NVT ensemble 
with the use of a Berendsen thermostat. The potentials 
are truncated at 9 A. The corrections to the long range 
Coulombic interactions are taken into account with the 
Ewald method. The timestep used is 1 fs. 

The density of water p is calculated on assuming that 
the excluded volume due to the confining spheres is ap- 
proximately Vexci = ''^NmO'om/^^ ^^i^ ^^-"^ ^ given 
number of water molecules 

- Vexcl Wmol 



where L is the boxlenght, Na is the Avogadro number 
and Wmol is the molar volume. In the starting config- 
uration the water molecules are arranged with the cen- 
ter of mass in a cubic lattice configuration where there 
are voids corresponding to the random positions of the 
matrix spheres. In order to have enough space for the 
spheres the boxlength is fixed to L = 19.25 A. Then for 
a given density we place in the box a number of water 
molecules according to Eq. |4l The range of temperatures 
spanned is 200 K < T < 350 K and the density range is 
0.75 g/cm? < p < 1.02 gjcrv? . For each density the sys- 
tem is melted at T = 500 K and equilibrated at several 
temperatures in the range indicated above. The longest 
runs, corresponding to the lowest temperatures investi- 
gated, lasted 10 ns. The soft spheres positions are the 
same for all the states simulated. We checked for an iso- 
chore that changing the disordered configuration of the 
soft spheres did not significatively alter the pressure val- 
ues (see fig. [TU)) . 

In Fig. [T] we show a snapshot of an equilibrated con- 
figuration for T = 300 K and p = 0.85 gjcrn?. It is 




FIG. 1: Snapshot of an equilibrated configuration of the sys- 
tem at T = 300 K and p — 0.85 gjcw? . Black spheres: 
hydrophobic solute; gray spheres: oxygen atoms; lightgray 
spheres: hydrogen atoms. 

evident the depletion regions around each hydrophobic 
sphere and the uniform distribution of water molecules 
elsewhere. 



III. STRUCTURAL PROPERTIES 

The site-site pair correlation functions have been cal- 
culated averaging the equilibrated configurations to see 
the effect of changing the density of water. The results 
for the water-water functions at ambient temperature are 
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reported in Fig. [2]|4]for different densities and compared 
with the bulk for p — 1 g/cnv'. The trend is similar to 
the one found in aqueous solutions of rare gase o^^i^°'^^ , 
considering that the decrease of the density in our case 
is equivalent to an increase of the solute concentration. 
The raise of OO peak in Fig. [2] with decreasing density 
is the signature of the enhancement of the structuring 
effect when less molecules are present and they can more 
easily approach one another. This is also evident in the 
behaviour of the goHif) in Fig.[3]and gHH{f) in Fig- 111 
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FIG. 2: Site-site pair correlation functions goo{r) at T = 
300 K for the confined system at densities p — 1.0 (dotted 
line), 0.90 (dashed line), 0.75 (long dashed line) g/cm^ com- 
pared with the bulk at T = 300 K and p — 1.0 g/crrC' (solid 
line) . 




FIG. 3: Site-site pair correlation functions goHij) at the same 
temperature and densities as Fig. [2] 

On the contrary on decreasing density the cages formed 
by water around the confining spheres become less struc- 




FIG. 4: Site-site pair correlation functions gHH{r) at the same 
temperature and densities as Fig. [2] 



tured as can be seen in Fig. [5] where the water-sphere 
radial distribution functions are reported for two differ- 
ent temperatures and different densities. The first peak 
of the goAiif) becomes less pronounced and shifts to- 
ward higher distances, a signature of an increase of the 
hydrophobic repulsion at lower density. In the upper in- 
set of Fig. O we show the number of nearest neighbors 
noM defined as follows 



noM = 47rpco 



goM{r)r^dr 



(5) 



where rmin is the value of the interatomic distance at 
which the first minimum in goAiir) is located, co is the 
concentration of oxygen atoms and p is the total density. 
We see that number of nearest neighbors noM reported 
in the top portion of the figure goes down as the water 
density decreases. In Fig. [5] it is also shown the effect 
of the decreasing temperature. The cage becomes more 
ordered and more water molecules are present on the av- 
erage in the first shell around the soft spheres at high 
density. On lowering the density the uom reaches an 
asymptotic value which appears to be independent from 
the temperature. 

The behaviour of gHAiir) also reported in Fig. [5] shows 
a similar trend as function of temperature and density. 
The first peak of gHM{r) is at the same position of the 
goKiif) one, this indicates that the soft sphere is located 
interstitially in the hydrogen bond network equidistant 
on the average from the oxygens and hydrogens. We no- 
tice in gHM{f) at lower temperature for the higher den- 
sity the appearance of a shoulder after the first peak, indi- 
cating that some of the 0-H bonds point radially toward 
the second shell. Looking at the second peak of f;oM('') 
and gHM{T) we notice that the water molecules in the 
second shell are preferentially oriented with the oxygen 
atoms towards the sphere and the hydrogen atoms to- 
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wards the water similar to the trend found in aqueous 
solution of water and rare gases^^iiSiil. 
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FIG. 5: Pair correlation functions goAiir) (middle panel), 
gHM(r) (bottom panel) for T = 300 K on the left and 
T = 220 K on the right and densities p = 0.90 (solid 
line), 0.85 (dotted line), 0.80 (dashed line), 0.75 (long dashed 
line) g/crrC' . In the top panel number of oxygen nearest neigh- 
bors noM as function of density for T = 300 K (full squares) 
and r = 220 K (filled triangles). 

In Fig. [S] to show the effect of the temperature on the 
structure of water we report the goo{f) for T — 300 K 
and T = 220 K and p = 0.90 g/cm^. The radial function 
shows maxima and minima at the same position but more 
pronounced at lower temperature. We notice that the 
gooi"!^) for p = 0.90 gjcw? in the confined system is 
very similar to the goo^r) for p = 1.0 g/cm? in the bulk 
apart from a difference in the height of the first peak. 
This similarity is preserved at decreasing temperatures. 

To complete the analysis of the change of density on 
structure of water we consider the changes in the H-bond 
network. We define the hydrogen bond by geometric con- 
ditions. Two molecules are hydrogen bonded if the in- 
termolecular O ■ ■ O distance is less than 3.35 A and the 
angle between the intramolecular O — H vector and the 
intermolecular O ■ ■ O vector is less than 30°. 

In Fig. [7] it is reported the distribution function of the 
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FIG. 6: Site-site pair correlation functions gooij') for the 
confined system at density p — 0.90 g/cm^ and temperatures 
T = 300 K (long dashed hne) and T = 220 K (dotted line) 
compared with the bulk at density p = l.Q g/cm^ and tem- 
peratures T = 300 K (solid line) and T = 220 K (dashed 
line). 



average number of H-bond per molecule at two different 
temperatures and three densities. We have computed the 
average HB number for each istantaneous configuration 
and then we have calculated the distribution function of 
these average values over many configurations. We now 
discuss the results for the bulk at normal density com- 
pared with p = 0.90 gjcrr? in the confined system since 
they appear to have a similar structure as shown above. 
It is evident that decreasing water density in confinement 
has the effect of reducing the average number of hydro- 
gen bonds in comparison with the bulk even when goo {t) 
are similar. At ambient temperature in confinement the 
distribution becomes broader and its peak position shifts 
to lower value with respect to bulk. These effects are en- 
hanced by decreasing the density. At lower temperatures 
the H-bond distributions becomes sharper while the peak 
position moves toward a value close to 4 

The confinement however preserves the geometry of 
the H-bond network. This can be deduced from the dis- 
tribution of the angle Q formed by the vectors joining 
the oxygen atom of a central molecule and the oxygens 
of two nearest neighbors water molecules which are H- 
bonded. This distribution is reported in Fig [51 It is 
evident that the trend is very similar between the bulk 
and the confined case. On lowering the temperatures the 
H-bond network approaches more closely the tetrahedral 
ordering with a sharper distribution around the angle 
104° of the TIP4P potential, while the number of inter- 
stitial molecules decreases. We note that the confined 
p = 0.90 g/cm? system has a distribution similar to the 
bulk p ~ 1.0 g/cm? system. 

In analogy with previous work on water confined be- 
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FIG. 7: Distribution functions of the average number of hy- 
drogen bond per molecule at temperatures T = 300 K and 
T = 220 K for the bulk at density p = 1.0 gjcrr? (solid line) 
and for the confined system at densities p = 0.90 gjcw? (long 
dashed line) and p — 0.75 g/cm^ (dashed line). 
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FIG. 8: Distribution functions of cos(ff) (see text for defini- 
tion) at temperatures T = 300 K and T = 220 K for the bulk 
at density p = 1.0 gjcvn? (solid line) and for the confined 
system at densities p — 0.90 gjcw? (long dashed line) and 
p = 0.75 gjcrr? (dashed line). 



tween hydrophobic plane walls^^'^^ we find a reduction of 
the average number of H-bond. In our case the H-bond 
network does not appear to be deformed by the presence 
of the matrix. In the hydrophobic shells around the soft 
spheres water molecules maintain the H-bond network, 
but since they cannot form H-bonds with the neutral soft 
spheres the average number of hydrogen bonds is smaller 
compared with the bulk. As a general trend we observe 
that confined water seems to be equivalent to bulk water 
at an higher density (see Fig. [HI), it is also evident that 



there are not dramatic changes in the network apart from 
a decrease of the hydrogen bond average number. 



IV. CALCULATION OF THE SPINODAL LINE 

As we have shown in the previous section the soft 
sphere matrix does not induce large changes in the struc- 
ture of water. We investigate now how the thermody- 
namical behaviour of confined water is changed with re- 
spect to the bulk and we will see that for thermodynamics 
important changes take place. 

To study how the confinement in the hydrophobic ma- 
trix affects the thermodynamic stability of liquid water 
we calculate the liquid branch of the isotherms of our 
system upon cooling. From the isotherms the limit of 
mechanical stability can be determined in the framework 
of mean field theories from the divergence of the isother- 
mal compressibility. The singularity points where 



dP 
'd'p 



= 



(6) 



define the spinodal line. They can be obtained search- 
ing for the minima of the Pt{p) curve. We consider 
here the liquid spinodal. In bulk water the liquid spin- 
odal^ii^i^iiaiiiii^ii^ has been found to extend to the re- 
gion of negative pressures. 
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FIG. 9: Isotherms for densities ranging from p = 0.75 g/cm^ 
to p — 1.0 g/cm?' for temperatures ranging from T = 350 K 
to T = 200 K. 

In Fig. [5] we show the isotherms of the system obtained 
upon cooling for densities ranging from p = 0.75 g/cm? 
to p — 1.0 g/cm? and temperatures from T = 350 K to 
T — 200 K. Below T — 280 K minima start appearing in 
the curves in correspondence with p = 0.8 g/cm? . In the 
bulk^ the minima appear between p = 0.85 — 0.87 g/cm^. 
From Eq. [5] we see that the minima of the isotherms 
represent the limit of stability of the liquid. The cor- 
responding Ps{T) spinodal line is plotted in FigfTUl to- 
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gether with the isochores of the system. In Fig[TO] we 
report for p = 0.95 g/cm?' isochore also the resuhs ob- 
tained for a different configuration of soft spheres. As 
it can be seen changing the soft spheres position does 
not alter significatively the values of the pressure. From 
the isochores also the behaviour of the TMD line can be 
extracted. In fact along the TMD line the coefhcient of 
thermal expansion 



1 fdV 



goes to 0. 

The thermal pressure coefficient 



(7) 



iv 



dP 
df 



(8) 



is connected to the coefficient of thermal expansion by 
the following equation 
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FIG. 10: Pp{T) of the confined system isochores for several 
values of p in g/cm? , spinodal line (long dashed with points) 
and TMD curve (solid line) . * Open circles represent the p = 
0.95 g/cm? isochore obtained from a diflterent configuration 
of the matrix of soft spheres. 

Therefore the TMD points lie on the line connecting 
the minima of the isochores. In the Fig. [10] also the TMD 
line is drawn. We observe, down to the lowest temper- 
ature investigated a non retracing spinodal and corre- 
spondingly a TMD that bends on approaching the spin- 
odal. This behaviour is similar to both non polarizable^ 
and polarizablei^ bulk TIP4P water. The interesting fea- 
ture is represented by the shift induced by confinement on 
both the spinodal and the TMD with respect to the bulk. 
To give a more quantitative insight on this displacement 
we plotted in FiglTT]the isochores and the TMD for the 
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FIG. 11: Pp{T) of the bulk system isochores for several values 
of p in g/cm^, TMD curve (long dashed line). The isochore 
of the lowest density is close to the spinodal. 
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FIG. 12: TMD and spinodal of the confined system (bold 
lines with filled circles) and TMD and isochore of the lowest 
density for the bulk system (bold lines with filled diamonds). 



bulk TIP4P. Some of the isochores have been calculated 
here as they are not available in literature. By comparing 
the spinodal and TMD curves of the bulk and confined 
system, FigfT2l we note that both the spinodal line and 
the TMD for water under hydrophobic confinement are 
shifted toward higher pressures. The TMD line shifts also 
to lower temperatures and correspondingly its curvature 
broadens. In a different hydrophobic confinement^Si, as 
already quoted in the introduction, a negative shift of 
40 K has been observed with respect to the bulk for the 
TMD, but no shift in pressure nor changes in the curva- 
ture were found. Interestingly our shift in temperature 
is also roughly 40 K. 



7 



V. CONCLUSIONS 

We performed a MD study of TIP4P water confined 
in a rigid matrix of soft spheres. The geometry of the 
hydrophobic environment in our case is different with 
respect to the previous studies mentioned above where 
water is confined between plane surfaces. We observed 
that the effect on water structure is similar to the case 
of a solution of small apolar solutes when it is taken into 
account that the decrease of water density is equivalent to 
an increase of solute concentration. Moreover we found 
that volume excluded effects reduce the average number 
of hydrogen bond, The hydrogen bond network however 



is preserved at variance with the case of water confined 
in hydrophobic plates. 

In spite of the preservation of the network of water 
and the not large changes in its structure we observe 
an important shift of the limit of stability both in pres- 
sures and temperatures and also a different shape in the 
shifted TMD with respect to the bulk. The tempera- 
ture shift similar to the one observed by Kumar et ali^. 
These simulation indicate only a weak dependence of wa- 
ter properties on the confining potential^''. It would be 
therefore valuable to understand to what extent an anal- 
ogy can be drawn. 
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